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ABSTRACT 



The force-free limit of magnetohydrodynamics (MHD) is often a reasonable ap- 
proximation to model black hole and neutron star magnetospheres. We describe a 
general relativistic force-free (GRFFE) formulation that allows general relativistic 
magnetohydrodynamic (GRMHD) codes to directly evolve the GRFFE equations of 
motion. Established, accurate, and well-tested conservative GRMHD codes can simply 
add a new inversion piece of code to their existing code, while continuing to use all the 
already-developed facilities present in their GRMHD code. We show how to enforce 
the E • B = constraint and energy conservation, and we introduce a simplified general 
model of the dissipation of the electric field to enforce the B 2 — E 2 > constraint. We 
also introduce a simplified yet general method to resolve current sheets, without much 
reconnection, over many dynamical times. This formulation is incorporated into an 
existing GRMHD code (HARM), which is demonstrated to give accurate and robust 
GRFFE results for Minkowski and black hole space-times. 

Key words: black hole physics, galaxies: jets, gamma rays: bursts, MHD, stars: 
neutron, magnetic field, outflows, methods: numerical 



1 INTRODUCTION 

General relativistic force- free electrodynamics (GRFFE) 
is the low-inertia limit of general relati vistic mag- 
netohy d rodyna mics (GRMHD) (see, e.g., iKomissarovl 
' l2002bl. l2004aTl . Neutron star and black hole mag- 
netospheres exhibit regions of space that are very 
nearly force-free, and self-consistent moderately realis- 
tic quasi-analytic solutions exist that describe the ideal 
MH D or force-free envi r onme n t of such systems (see , 
e.g.. iBlandford fc Znaiekl Il977t IContopoulos et alJ Il999t 
[ Goodwin et alJ 12004 iBeskin fc Nokhrinal 120051: iGruzinovl 
2005). It is generally difficult to solve the GRFFE equa- 
tions to find even stationary solutions, except with sim- 
plified assumptions, that apply to astrophysical systems. 
For example, despite recent progress in studies of neu- 
tron star magnetospheres, no self-consistent analytic solu- 
ti on considers general relativ istic effects. Also, the solution 
of lBlandford fc Znaiekl <ll977f) is the only self-consistent an- 
alytic force-free solution that has a realistic Poynting jet. 
Also, analytic solutions typically assume stationarity and 
axisymmetry and so rarely address the global or local sta- 
bility of the solutions against time-dependent or nonaxisym- 
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metric modes or stability against reconnection in current 
sheets if present. 

One simple technique to seek stationary solutions and to 
study time-dependent stability is to directly numerically in- 
tegrate the GRFFE equations of motion. Both the GRMHD 
and GRFFE equations of motion can be written as a set 
of conservation laws that can be directly integrated. Con- 
servative numerical G RMHD methods, such as of HARM 
jGammie et alj l2003a1 , use so-called "primitive" quantities 
(P) to define so-called "conserved" quantities (U), fluxes 
(F), and source (S) terms. The temporal integration is deter- 
mined by solving the set of "conservation" equations, which 
can be written as 

d^ + b (P) ' (1) 

where p labels the conservation equation and i is the spa- 
tial index. The set of source terms (S p ) accounts for the 
connection coefficients or other sources of mass-energy- 
momentum. Thus conservation is explicitly true as long as 
the source terms vanish. For an axisymmetric, stationary 
metric and as written in HARM, energy and angular mo- 
mentum are explicitly conserved to machine error because 
the source terms vanish. HARM has been successfully used 
to study black hole accretion flows, winds, and Poynting 
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jets jGammie et alJl2003at Ig ammie. Shapiro, fc McKinnevI 
l2004lMcKinnev fc Gammidl2004UMcKinnevll2005allbld) . 

While GRFFE is mathematically just the no-inertia 
limit of GRMHD, numerical truncation errors limit the use 
of GRMHD codes in evolving systems with regions where 
the magnetic ene rgy density greatly exceeds the rest-mass 
energy (see, e.g ., iMcKinnev fc Gammid Eooit iKomissarovl 
l2004bl. l2005allbK Conversely, the GRFFE equations of mo- 
tion do not describe the rest-mass motion along field lines 
or the thermal energy of the particles. So long as the iner- 
tia of particles remains negligible, the GRFFE equations of 
motion properly describe the magnetic field geometry and 
motion. 

A GRFFE code can be used to study neutron star and 
black hole magnetospheres. For example, the origin of the 
collimation and stability of astrophysical jets remains un- 
explained. A GRFFE code can be used to examine the ori- 
gin of the collimation and the stability of any model for a 
Poynting-dominated jet. 

Existing jet solutions may not be sta ble and self- 
collimation of isolated je ts may not work iSpruitl Il996t 
IOkamotoll 999 . 200rj. 12003) . Collimation by a predominately 
toroidal field (hoop stress) may lead to a pinch/kink in- 
stability, whereas poloidal collimation works only if the 
jet is surrounded by an extended dis k wind dSpruidri996t 
iBeeelman fc Lilll994l ; rBegelmanlll998l) . The solution to the 
collimation problem may be that the disk wind collimates 
the Poynting jet. By using both GRMHD and GRFFE nu- 
merical models, the origin of collimation, acceleration, and 
stability of jets can be examined. Thus it is more efficient 
to have a single code be able to perform both GRMHD and 
GRFFE studies. 

Even a black hole system a ccreting a thick disk 
has a force-free magnetosphere iMcKinnev fc Gammid 
120041) . Such a ma gnetosphere quantita ti vely a grees with 
the solution by iBlandford fc Znaiekl (Il977t) for low 
black hole spins and agrees qualitatively for hi g h spin s 
iMcKinnev fc Gammid 12004 IKomissarovl l2004bl. l2005af) . 
Gamma-ray bursts (GRBs), active galactic nuclei (AGN), 
and x-ray binaries probably exhibit both a Poynting- 
dominated jet and a jet fr om the disk, as suggested by 
GRMHD numerical models iMcKinnevll2005bT lcll . In partic- 
ular, of all jet mechanisms, the Blandford-Znajek driven jet 
is the o nly one that can clearly produce an ultrarelativis- 
tic jet jMcKinnevll2005ari . Thus is it important to study 
the Blandford-Znajek process in detail. A GRFFE code can 
help determine the stability of the Poynting-dominated jet 
generated by the Blandford-Znajek effect. For example, a 
GRFFE code was used to study the pure monopole version 
of t he Blandford-Znajek solut ion, which is found to be sta- 
ble (IKomissarovl 1200 ll l2002ah . However, general Poynting- 
dominated jets produced by the Blandford-Znajek effect 
may be viol ently unstable to pinch and kink instabil ities 
(see, e.g.. IlI EoOoI but see also iTomimatsu et al.|l200ll) . A 
GRFFE code proves valuable in stability analyses by avoid- 
ing the difficul ty of analytically crossing the light cylinder 
(see, e.g.. ITomimatsu et al.ll200ll) . which manifests itself as 
a singularity in the Grad-Shafranov equation for stationary 
solutions. 

The goal of this paper is to formulate the GRFFE equa- 
tions of motion such that a conservative GRMHD code can 
be used with small modifications. Present numerical meth- 



ods directly evolve the electric and magnetic fields (primi- 
tive quantities are {E, B}, the field-evolution approach) and 
have to explicitly check that the velocity is less than the 
speed of light (i.e. magnetic energy density is greater than 
electric energy density, B 2 — E 2 > 0) since the electric and 
magnet ic field are evolved with no constraint on thei r evo- 
lution jKomissarovll2002bl: iKrasnopolskv et al.1 12005). The 
method we describe evolves the drift velocity and magnetic 
field (primitive quantities are {v, B}, (velocity-evolution ap- 
proach), explicitly breaks down when that constraint is vio- 
lated and this ensures a physical and causal evolution. 

We show that this velocity-evolution approach can 
be directly used by GRMHD codes, and explicitly dis- 
cuss how this is done for conservative numerical meth- 
ods. Many GRMHD codes have recent l y been devel- 
oped jKoide et all Il999l: I K^mjssarov| ^^99^jGammi^^^alJ 
| 2003al: iDe Villiers fc Hawleyl 120031 IShibata fc Sekiguch: 
120051: lAnninos et alJ 120051 : lAnton et alJ 120051) . Our dis- 
cussion applies to any GRMHD method, but focuses on 
conservative-based codes. We show that a separate GRFFE 
code does not have to be developed, and many of the tools 
and methods used to make the GRMHD code perform well 
carry directly over to the GRFFE version. 

§|2]shows how a conservative GRMHD code can be used 
to evolve the GRFFE equations of motion. Simplified models 
of dissipation in GRFFE are discussed in order to handle 
regions here the electric field dominates the magnetic field. 

§ includes a series of standard tests of the GRFFE 
code. 

§0discusses models of black hole magnetospheres. This 
section demonstrates the usefulness of the GRFFE equations 
of motion and the model of dissipation. 

§ |5] summarizes the results of the paper. 

The GRMHD notation follows iMisner et alJ Jl973t) and 
we use Heaviside-Lorentz units unless otherwise specified, 
which is like Gaussian units without the 4n. For exam- 
ple, the 4-velocity components are u M (contravariant) or 
Uft (covariant). For a black hole with angular momentum 
J = jGM 2 /c, j = a/M is the dimensionless Kerr parameter 
with — 1 < j < 1. The contravariant metric components are 
<7 M " and covariant components are g M „ . The comoving energy 
density is b 2 /2 , where is the comoving magnetic fi eld. See 
iGammie etTafl d2003al) : lMcKinnev fc Gammid' (120041) for de- 
tails. 



2 FORCE-FREE ELECTRODYNAMICS 

This section shows how a conservative-based algorithm de- 
signed for solving the ideal GRMHD equations of motion 
based upon the velocity and a magnetic field can be used 
to solve the GRFFE equations of motion. In particular, the 
electric and magnetic field are used to derive the velocity 
of the frame in which the electric field vanishes. This al- 
lows a conservative ideal GRMHD numerical code to evolve 
the force-free equations of motion by simply modifying the 
inversion from conserved quantities to primitive quantities. 

The force-free electrodynamics equations of motion are 
the 8 Maxwell equations: 

V MJ F M " = -J", (2) 
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where F' 1 " are the components of the Faraday tensor and 
J M are the components of the current density, and 

av 



0. 



(3) 
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where F 

Faraday (or Maxwell) tensor, F^ v 
{— l/y/— g)[a/3"f8], and e a ps-y = y/—g[a/3'y8], where brack- 
ets denote the unit length completely antisymmetric tensor. 
These 8 equations define 6 evolution equations and 2 differ- 
ential constraints. 

By the antisymmetry and duality of the Faraday and 
Maxwell tensors, they be expressed in component form as 



a0 



F 
and 



frfE" 



hB^rise 



p p = - hr) <* B e + hrfB a 



(4) 



(•») 



where / and h are arbitrary independent constants 
that we set t o be / = h = 1 consistent with the 
conventions in |Mjsngr_et_alJ <ll973l) . and see also, e.g., 
iBaumgarte fc Shapirol tOOJ) . Here E a = r ! 0F c " 3 , B a = 

* 0a 

rjg F , and ry M is an arbitrary 4-vector. This gives that 
F^F^ = 4 J E fl B fl (-r ? 2 ) and F^F^ = 2{B 2 - E 2 ){-~r) 2 ). 
For a time- like jf" , the sign conventions are as in special 
relativity. 

The electromagnetic stress energy tensor is constructed 
from the Faraday tensor and must be quadratic in the field 
strengths, symmetric, and is divergenceless in vacuum by 
Maxwell's equations. This unique tensor is 

1 



(6) 



where g^ v are the components of the metric. The duality be- 
tween the Faraday and Maxwell tensors and the definitions 
of and give that 

* (77V + -PH (7) 



-rf) (B"B V + E»E V ) 



rj a Ef}B K 



U lfCt0K . V liOL0K 

Tj € + Tj € 



where the projection operator P M " — (— r/ 2 )g' J ' v + rfrf '. 
For example, the energy density in frame rf is r\^r\ v T tlv = 
r] 4 (E 2 + B 2 )/2, which for a time-like rf is the same ex- 
pression as in special relativity. If the electromagnetic field 
is the only source of stress-energy, then equations (5) are 
equivalent to the energy-momentum conservation equations 



V M T M " = = 0, 



(8) 



where only 2 of the 4 components of equation JHJ are inde- 
pendent because equations (JHJ implies 



f F V" =4 J B fl B M (-r ? 2 ) = 0. 

For more details see iKomissarovl d2002bll2004ah . 



(9) 



2.1 Inversion 

Conservative numerical GRMHD methods, such as of 
HARM, operate primarily on so-called "primitive" quan- 
tities (P): fluid density, fluid internal energy, coordinate 
fluid velocity, and the lab-frame coordinate magnetic field. 



HARM uses the primitive quantities to define so-called "con- 
served" quantities (U) and the fluxes (F), which are both 
closed-form operations. These conserved quantities can be 
evolved forward in time, but then an inversion to primitive 
quantities is required to easily define the fluxes and other 
required quantities to determine the next update. 

No closed form solution appears to exist in GRMHD, so 
most inversion methods rely on iterative procedures such as 
Newton's method. This approach can be used in force-free 
electrodynamics, but may not always work due to known 
problems (such as poorly conditioned or singular Jacobians) 
with Newton's method. However, a closed-form solution does 
exist for force-free electrodynamics, as described below. 

The only conserved quantities of relevance in force-free 
electrodynamics are the lab-frame momentums (the en- 
ergy evolution equation giving the /i = t term is actually 
redundant) and the lab-frame magnetic field _B M (only 3 
components are independent). If one could obtain E^ , then 
one could reconstruct the Faraday from equation @ or any 
other quantities from E^ and B M . It is straight-forward to 
show that if E a B a = 0, then 



E a = e 



BpT^ Vs /(B 2 V 4 



(10) 



as shown by a substitution of equation @ and equation (g( 
into this expression. Only the spatial components of J^T^ 
and _B M are needed if one chooses a special form of rf" . Notice 
that E a B a — is explicitly true, therefore the degeneracy 
condition of E ■ B = can be preserved to machine error 
regardless of the truncation error in T or B. Also notice that 
equation 11011 projects out a component perpendicular to 
time, the spatial field, and the momentums. For a fixed field 
B l , only 2 components of the electric field are independent. 

One may choose to have rf = — 1 such that rf" only de- 
pends on the metric. One choice is = {—a, 0, 0, 0}, where 
a = l/v 73 ^"- Then rf = (l/a){l, -/?}, where j? = a 2 g u . 
This defines a zero angular momentum (ZAMO) frame. This 
choice of -rf" makes it possible to only require T\ and B l to 
obtain E a in equation (1 1 ( )fi . 

Another interesting choice for is to have E" = 
r\y.F^ v — 0. In the ideal GRMHD equations of motion, for 
a fluid velocity tt M , this choice corresponds to the electric 
field in the comoving frame w M being e" = u^F^" — 0. Then 
rf = — 1 since the fluid velocity is time-like. In force- free 
electrodynamics, there is no unique 4-velocity that satisfies 
e v — 0, but one such frame is constructed below that is 
uniquely always time-like with a minimum Lorentz factor 
with respect to the frame with 4-velocity rf . 

As shown next, any two frames with 4- velocities u M and 
?7 M can be easily related to determine a 4-velocity of the 
frame in which e v — for the Faraday defined in terms 
of ?7 M . Thus the metric and Faraday alone determine a 4- 
velocity that allows one to use the ideal GRMHD Faraday, 



F 



a0 



a{3~f& 



and Maxwell, 



a 



a i0 i pii 

u b + u b 



(11) 



(12) 



where b" = lt M F . This formulation and sign conventions 
are the same as used in HARM. In HARM, a new 4-velocity 
is introduced that is unique by being related to a physical 
observer for any space-time and has well-behaved interpo- 
lated values, 
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(13) 



where 7 = —u a ri a and so u = t/ck- This additional term 
represents the spatial drift of the ZAMO frame defined ear- 
lier. One can show that 7 = (1 + g 2 ) 1 ^ 2 with q 2 
To obtain the 4- velocity, notice that 

= e a = -{vpE ) n a + £ Q + v B im e a ^ s , 



gijU % v? . 



(14) 

where = —up/{u a r]a.). The most general form of the 4- 
velocity that satisfies the above is 



Vf} = G 



+K 



-r] 2 B 2 



H 



10 



^Fe 2 



+ L 



Be 



(15) 



where G, H, K, and L are functions to be determined. 
Each term represents one of four orthogonal directions, and 
when multiplied by arbitrary functions represent the most 
general solution. The only nontrivial term is the antisym- 
metric product between the last term in equation I14H and 
the first term in equation 1151 . Substitution of this v 13 
into equation Q140 gives that G = 1. Since vprj 13 — —1, 
then H = 1. With up — ■yvp, then u 2 — — 1 gives that 
7 = ^rfB^JJJl - K 2 )B 2 - E 2 ), which simply defines 
7 = — u a r) a as the Lorentz factor between the two frames. 
All terms proportional to K are orthogonal to each other 
and to E a , and so in general K — 0. Hence, u a E a — 0. To 
determine L, notice that b a in F in equation 112H can be 
written as 



P$B fl 
7 



(16) 



where Pp = dp +u a up is the projection tensor. Thus b a has 
terms proportional to B a and u a . So extra terms added to 
u a proportional to B a vanish due to the antisymmetry of 
the Maxwell and so do not contribute to the stress-energy 
tensor or the equations of motion. Thus, the function L pa- 
rameterizes the arbitrary velocity component along a field 
line, and we choose L = to minimize the Lorentz factor of 
the frame. Hence, u a B a = and so b a — B a /j. With this 
choice of L, if B 2 — E 2 > 0, then the frame with u M is always 
time-like. Thus the frame defined by the 4-velocity 



B 2 



B 2 



E 2 



r)pE 1 Bs 



B 2 



(17) 



is time-like for any force-free electrodynamic solution. Thus, 
by construction, we have shown that in force-free electrody- 
namics that it is possible to boost into a time-like frame 
where the electric field vanishes and thus the Poynting flux 
vanishes (see also, e.g. lKomissarovll2002br) . This also shows 
that force-free electrodynamics is a causal limit of GRMHD 
as long as B 2 — E 2 > 0. This 4-velocity also represents a 
unique covariant definition of the "field-line velocity," and 
also describes the field- line velocity even in ideal MHD. 

A GRMHD code may only need the coordinate lab- 
frame 3- velocity. Since it* = 777*, then for the earlier defined 
ZAMO frame , the coordinate lab- frame 3- velocity is given 
by 



: [ijk]EjB k 



v = — = -$ +a 

u y—gB 2 



(18) 



The second term in equation 1171 represents the purely spa- 
tial "E x B drift." Note that there is no evolutionary con- 
straint on T' 1 " that forces B 2 - E 2 > 0, and when this 
is violated the force-free model is no longer physical. The 
value of tr coincides with the "field rotation frequency" 
flf = F t r/F r< f, = Fto/Fgip for stationary axisymmetric flows 
for which v r = v° — 0. 

Now the inversion from {T^B^ ~* {E^B*} and then 
— > {v t ,B 1 } completely defines equations 111! and 1121 for a 
general relativistic force-free electrodynamics evolution us- 
ing a conservative GRMHD code. 

Obviously this formulation can be also used to study 
special relativistic models as well. Notice that in special rel- 
ativity that the derived velocity expression reduces to the so- 
called "electromagnetic 3- velocity" v = E x B/B 2 = S/B 2 , 
where S = E x B is the Poynting flux. 

As used in HARM, this formulation preserves V ■ B = 
and E ■ B = to round-off error for both the GRMHD and 
GRFFE equations of motion. Notice that this differs from 
other formulations t hat only preserve V B = and E-B = 
to truncation error f Komissarov 2002bl. I2004H) 

The fact that the field-evolution approach does not 
break down might be considered an advantage when seek- 
ing stationary solutions. In such a case the evolution may 
have regions that only transiently have B 2 — E 2 < and 
the integration can pass smoothly through this region into 
a physical solution space. A related advantage of the field- 
evolution approach is that Runge-Kutta temporal evolution 
can recover the correct temporal order of accuracy. That 
is, without characteristic interpolation, HARM uses Runge- 
Kutta to time step to achieve higher order temporal accu- 
racy. Runge-Kutta is only first order accurate for the first 
substep, but after all substeps are completed, the method 
is accurate to arbitrary order. The velocity-evolution ap- 
proach can yield an unphysical result for a substep and be 
unable to continue or treat the result as a violation of force- 
free electrodynamics, while the field-evolution approach can 
avoid such first order errors and recover to arbitrary order 
accuracy. However, we have found the velocity-evolution ap- 
proach to be sufficient. This issue of falling outside the light 
cone is the same issue one encounters when evolving the 
GRMHD equations of motion, and in that case the author 
knows of no method that does not require the velocity at 
some point during the integration, so the issue is treated as 
a generic one that simply requires a more accurate integra- 
tion. 



2.2 Currents 

This section shows that the currents can be computed with- 
out time derivatives, which is numerically convenient to 
avoid storing data at previous times. In ideal MHD or force- 
free electrodynamics there are many dependent ways of 
equally describing the same physics. Researchers often in- 
voke, such as in discussions of current closure, the current 
and the magnetic field as a more intuitive set of quantities 
than the electromagnetic fields. The current could be com- 
puted from 



J° 



(19) 



which apparently requires time derivatives. In force- free elec- 
trodynamics, however, the current (like the 4-velocity) sits 
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in the null space of F al3 , i.e. J a F^ = 0. Hence, J a E a = 0. 
The same analysis as in section |2~T1 leads to the same result 
for J a except the function L is no longer arbitrary, where 
here 



J 3 



(20) 



where p qtV = J a ri a , which only actually requires spatial 
derivatives if 7)^ = {—a, 0, 0, 0}. Plugging equations Q into 
equations 1191 and contracting Jp with B 13 gives 



J B P = [E a r,g, a - r) a E p , a - (B 7%e/3 



(21) 



where the first and last terms only require spatial derivatives 
with the chosen r\^. Using 



F;3 =0, 



(22) 



and contracting with E a , the second term in equation 1211 1 
can be written as 



B p r, a Ef3. :a = B a E p Vp . a - Ep. a E^ s e 



(23) 



Now, for the chosen r/^, each term in equation 12111 only 
requires spatial derivatives. So the current can be written, 
only actually requiring spatial derivatives, as 



= Pq.V [V 



e a ^ s V0 E y B 5 \ +Ba (J B B g 



B 2 ) \ B 2 

where finally 

J p Bp = E a B 13 (fi P , a + rj a ., B ) 

+ (B a . B 1 - Ec-pKy) Us^"' 5 



(24) 



(25) 



where the asymmetry between E a and B a just expresses 
the asymmetry in Maxwell's equations for J a . One can use 
equation 1171 to write 



J = U 



(26) 



That is, there is a perpendicular drift current (Jl) and a 
field-aligned current (J\\), i.e. 



J a = J? + Jii 



(27) 



Equation 1261 is Ohm's law in dissipationless GRFFE. In 
the special relativistic regime this reduces to a lab-frame 
current density of 

_ ExB(V-E) + B(B-(VxB)-E-(VxE)) ,„ Q . 
= b 2 ' ( ' 

where rj a J a = p q , v = —a J* = —V ■ E is the charge density 
in the frame moving with 4- velocity 77^. As in general, the 
special relativistic equation obviously has no time deriva- 
tives. 



2.3 Jump Conditions 

Conservative schemes are often based upon ID piece-wise 
constant Godunov methods that achieve higher than first or- 
der accuracy by relying on a one-dimensional interpolation 
from, e.g., grid cell centers to cell interfaces. The interpo- 
lated states are assumed to be approximately constant over 
the cell face, otherw ise a genera lized Riemann problem must 
be solved (see, e.g. lTorolll999l) . The ID Godunov scheme 



then explicitly treats the cell interface discontinuity appro- 
priately and reduces to a trivial form for smooth flows. How- 
ever, while the hydrodynamic equations of motion allow ar- 
bitrary left and right initial states, the electromagnetic field 
must obey general well-known jump conditions as a mani- 
festation of the Bianchi identities and the antisymmetry o f 
the Faraday tensor (see, e.g., chpt. 15 of lMisner et al.ll973f) . 
Alternatively stated, electrodynamics can be written as a set 
of Hamilton- Jacobi equations involving the vector potential 
and the electric field. Arbitrary interpolation of arbitrary 
quantities does not generally enforce these jump conditions. 
Thus, these jump conditions must be self-consistently en- 
forced by interpolating appropriate quantities in the appro- 
priate way. 

HARM directly operates on the magnetic field rather 
than the magnetic vector potential, so the scheme must 
enforce the electrodynamic jump conditions on the elec- 
tric and magnetic fields. For schemes that only interpolate 
in space and not time, one only requires the spatial jump 
conditions. The jump conditions across a one-dimensional 
space-like surface for a single lab-frame time (t) are ob- 
tained by integrating Maxwell's equations. For simplicity 
define ^—gE a = t 13 F a0 and y/—gB a = t 13 F 0a . Also, con- 
sider three arbitrary orthogonal space-like vectors A,B,C 
that describe a space-like volume and the time-like (often 
but not always Killing) vector t a = {1,0,0,0}. For the ho- 
mogeneous Maxwell equations, 



= 



, * lll> . 

-gF ),« 



(29) 



where >' 3 = [ijk]E k . First, let dE M = e f , a p- / t a dA f3 dB J be 
a 1-form 2-volume for a single lab- frame time (see, e.g., 
lLichnerowic3ll967t Il976l lAnildll989h . Then equation 
gives for the two spatial parallel components that 



-gE A d A \CAB) = 0, 



(30) 



where the C-direction is chosen as perpendicular to the sur- 
face and never is a sum implied for [CAB], which only 
gives the 3-signature for arbitrary, but fixed, A,B,C cor- 
responding to any 3 spatial directions. Choosing instead 
dE M = e l _ la0 ~ l dA a dB 13 dC 1 , one has for the contravariant field 
that 



Ac / V^B^dAdB = 0. (31) 

Likewise, for the inhomogeneous Maxwell equations, 
j - < 

where F %] 



-dE u 



[ijk]B k . With dE p 



-gB A dA = [CAB] 



(32) 

e fJ , a /3- y t a dA l3 dB'i, then 



-g.J B dAdC, 



(33) 



for a possible surface current J B = 5(C)K B , where upper 
(lower) A,B,C denotes the contravariant (covariant) com- 
ponents parallel to the surface. Choosing instead dE M = 
e fla/3l dA a dB ,3 dC J for space-like orthogonal A, B, and C, 
one has for the contravariant field that 



-gE°dAdB 



V^J^V^gdAdBdC, 



(34) 
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for a possible surface charge p QtV = r/^J^ = S(C)a q , v , where 
a q , v is the surface charge in the frame moving with 4- velocity 
rf" and the C-direction is taken to be across the surface. 

For an infinitesimal surface yJ~gdAdB, lengths yJ~gdA 
and y/—gdB, or for point values of the fields, these expres- 
sions reduce to 

A c Ea[CAB] = 0, (35) 
A C B C = 0, (36) 
and 

A c Ba = [CAB]K B , (37) 

A C E C = a w (38) 

which apart from the concern with y/—g, the contravari- 
ant vs. covariant components, and E, B vs. E, B is identi- 
cal to the Cartesian special relativistic expressions for the 
jump conditions. Notice that no distinction is made between 
bound and free charges or currents, so these expressions are 
generally true. 

The equations (I35H through 138H must be preserved at 
the 1-D cell interfaces. For centered schemes the continu- 
ities can be enforced by using the same interpolation stencil 
for the left/right interpolations to the cell interface from 
the left/right cell centers. This is the method chosen for the 
GRFFE version of HARM. These constraints on the fields 
form an implicit constraint on the drift velocity and the mag- 
netic field. Notice that equations 1371 and 13811 enforce no 
specific constraint unless there is an enforced surface charge 
or surface current. Also notice that these constraints must 
also be preserved in ideal GRMHD. 

Interpolating the electric and magnetic field can lead to 
unphysical interface states if B 2 — E 2 < is unconstrained at 
the interface. To avoid this, one can interpolate 7_E Q and 7 
separately and then reconstruct E a for a given interpolated 
7 at the interface. This guarantees that the interface state 
is physical and reasonably similar to the center states. 

An alternative scheme can be designed with a stag- 
gered grid that autom atically enforces these constraints 
( Del Zanna ct al .120031) . However, both methods involve the 
same number of interpolations since in their case they must 
interpolate the interface field components to the center be- 
fore performing the inversion from conserved to primitive 
quantities. They must also use a stencil that guarantees no 
discontinuities at their cell center. Thus both methods are 
equivalent. 

2.4 Energy Conservation 

The formulation above for the inversion in force-free elec- 
trodynamics only explicitly requires T* and not Tf , which 
is a similar feature to o ther methods tKomissarovl [2002b; 
iKrasnopolskv et all2005T) , Since the truncation error in each 
is independent, the conserved quantity associated with 
energy conservation (T*) is generally inconsistent with T*. 
Hence, energy is only conserved to truncation error. As for 
any numerical scheme that only conserves energy to trun- 
cation error, this error can be used to gauge the reliability 
of the integration. However, in steady-state problems this 
truncation error may be secular and build-up and lead to un- 
realistic solutions. It is fruitful to formulate the inversion to 



enforce energy conservation and compare the nonconserva- 
tive integration with a conservative one to gauge the actual 
effect of the truncation error. 

One requires a solution for T* as a function of an ar- 
bitrary set of three components of T^. Then an arbitrary 
choice can be made to move the truncation error from en- 
ergy to a momentum. For axisymmetric, stationary space- 
times the natural choice of momenta are the radial and 9 
momenta. This relationship between that only depends 
otherwise on B^ can be obtained from 

T^rf = r?^f~, (39) 
where equation (f 10pi also gives that 

Equation (M0> does not actually depend on T* for our choice 
of rf . Equation <l3"§t gives T\ in terms of T* and B l . Notice 
that contracting with only rf shows that rfTj} just reduces 
to i) a B !i T% = for v ^ t. This shows that of the three T/'s 
that only two are independent for a fixed B l . However, for 
an arbitrary magnetic field, the evolution of all three T*'s is 
required to avoid singular expressions. 

To replace any particular T/ with T t , one solves the 
quadratic equation 1391 for that spatial component in terms 
of T t and the remaining spatial components. This new set of 
effective T/'s can then be plugged into equation 1101 . Notice 
that in general the quadratic equation gives two solutions. 
This degeneracy is introduced by using the energy, which 
lacks directional information. In practice it is sufficient to 
use the solution closest to the one obtained originally from 
integration of only the spatial parts (T*). This procedure can 
be used to keep energy and angular mo mentum conserved to 
machine error, unlike in prior methods dKomissarovi r2002b; 
IKrasnopolskv et al.l2005T) . Care must be taken for numerical 
methods with a large truncation error in the energy, since 
the evolution of the spatial and temporal stress-energy com- 
ponents may be disparate. The origin of this larger trunca- 
tion error is the larger nonlinearity of the energy compared 
to the momentums. Lack of energy-momentum conservation 
can also be due to dissipative processes, as described in the 
next section. 



2.5 Dissipation in Force-Free Electrodynamics 

There is no evolutionary constraint in dissipationless force- 
free electrodynamics that preserves B 2 — E 2 > 0, whose 
violation is taken as evidence that the plasma being mod- 
elled would have a nonnegligible inertial back-reaction on 
the electric field. This typically occurs in current sheets or 
in regions where the inertia would restrict the field to have a 
velocity of 11 < c. A physical system responds by dissipating 
the electric field into other forms of energy. 

The dissipation of the electric field is determined by an 
Ohm's law. For magnetospheres with an ample supply of 
charges (i.e. not "charge-starved") the Ohm's law in force- 
free electrodynamics is well-approximated by a large conduc- 
tivity along the magnetic field and a v anishing conductiv- 
ity per pendicular to the magnetic field (iKomissaro vl l2004al. 
l2005bl) . This reduces Ohm's law to the condition E ■ B = 
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and the only perpendicular current is the drift current with 
velocity given by equation 11811. 

To model this dissipation, iKomissarovl j2005lJ) use the 
prescription that if B 2 — E 2 < 0, then they introduce a large 
cross-field conductivity. They also have to modify the drift 
velocity to keep v < c. The problem with this prescription is 
that there is no dissipation until v > c. Indeed, in currents 
sheets the rate of dissipation is related to the drift velocity 
and is allowed to be v — > c, leading to the fastest possible 
reco nnection ra te. 

iLvubarskvl J2005I) study the relativistic Sweet-Parker 
sheet and Petschek configuration and determine that the 
rate of reconnection is much less than th e speed of light , 
contrary to previous estim ates (see, e.g. iLvutikovl l2003t 
iLvutikov fc UzdenskvlEool . This suggests that dissipation 
should limit the drift of field into a current sheet. 

First, inertial losses due to the dissipation of the electric 
field in a relativistic flow, such as beyond the light cylinder 
of a rotating compact object, can be "captured" by limiting 
the Lorentz factor of the drift velocity. This is a similar 
approach taken in GRMHD numerical models in order to 
avoid significant numerical errors at large Lorentz factors. 
A simple prescription is to limit 7 = —u a n a such that 1 ^ 
7 2 ^ Imax- From — 1 = u a u ay this means replacing B 2 in 
equation 1181 (and the similar B 2 in the second term of 
equation ll7t with 



(-n 2 )E 2 B 2 



(41) 



only when 1 > -y 2 > "fm ax . This gives a limited 3- velocity of 



\ijk\EjBk 
—gP 2 



(42) 



that has a continuous transition between standard force- 
free electrodynamics and dissipative electrodynamics. The 
4- velocity given by equation 11711 is, by construction, limited 
to always be time- like and have 7 ^ 7 mm . The energy- 
momentum lost in such a limiting procedure is assumed to 
be gained by inertial mass in the form of, e.g., (rapidly lost) 
thermal or field-aligned kinetic energy that has no effect on 
the field. 

Second, current sheets dissipate due to advection of field 
into the sheet and subsequent reconnection and cancelation 
of the magnetic field. Even under exactly symmetric con- 
ditions, numerical round-off error quickly leads to random 
velocity drifts across the current sheet. This magnetic field 
advection can be limited or avoided by limiting the drift ve- 
locity perpendicular to the current sheet. For example, the 
3- velocity given by equation 14211 can be further modified to 
have a small or vanishing component perpendicular to the 
current sheet. That is, if n ! is the spatial normal vector of 
the current sheet plane at a particular time-slice, then we 
can set 



n v Qij — (J 



(43) 



within some infinitesimal region bounding the current sheet. 
This changes Ohm's law in equation 127H to have a vanishing 
conductivity perpendicular to field lines, i.e. the spatial part 
of J™ vanishes such that the "E x B drift" vanishes along the 
normal direction of the current sheet. This explicitly avoids 
significant reconnection by avoiding anomalous numerical 



drifts. This is useful in studying systems for which the effect 
of reconnection is uncertain and the current sheet may be 
stable. The above prescription can be generalized to set any 
arbitrary drift speed into the current sheet. Thus, physical 
models of the current sheet and reconnection speed can be 
included in force- free models as long as inertia plays no other 
role than in the sheet. 

At present this is only implemented for a priori known 
locations of the current sheet, although a general algorithm 
shoul d be similar to recon nection-capturing methods (see, 
e.g.. lStone fc Pringkjl200lh . In the tests below that have a 
current sheet, there are 4 numerical grid zones in a current 
sheet region that are forced to obey the above condition. 
This guarantees that the stencil, used by the reconstruction 
and other dissipative procedures in HARM, does not cou- 
ple quantities diffusively across the sheet. This also ensures 
that upon convergence testing that this numerical scheme 
to avoid reconnection in current sheets plays no role in the 
results. 

This formulation ensures that B 2 — E 2 > and that 
current sheets can be forced to be mostly stable, unlike in 



IKomissarovl l)2002bl ) : |Krasno_r2oIi jkv et all J200fih and is im- 
proved compared to lK^mis^arovMl2005bir 



2.6 Quasi-GRMHD and Stationary Force-Free 

While this paper describes a GRFFE formulation, since the 
same code also has a GRMHD formulation, one can imagine 
hybrid schemes that integrate the decoupled equations of 
motion in the stiff regime where b 2 /po ^> 1, where pa is 
the rest-mass density. All GRMHD numerical schemes have 
difficulties with this regime, so a decoupled evolution may 
prove useful for studying the first nonzero order effect of 
inertia in a force-free field. This hybrid method will be used 
in future work, but the method is outlined here. 

In the force- free limit, one may still retain the evolu- 
tion of the rest-mass and internal energy density with no 
back reaction onto t he field (see alsolMe s^el fc S hibatall994l : 
IContopoulol Il995l: IContopoulos et al.lll999h . The method 
is to evolve the full GRMHD equations of motion, but to 
determine the field-perpendicular velocity components and 
field from the force-free equations alone. In this case one can 
readily obtain the field-aligned velocity component from the 
GRMHD inversion. That is, in ideal GRMHD in general, the 
fluid velocity may be broken into a field-perpendicular and a 
field-aligned velocity or equivalently into an electromagnetic 
and a matter velocity, 



a , a a. .a 

l ± + U|| = U EM + U MA . 



(44) 



Also of interest is the fluid motion in a stationary force- 
free field, which allow s a study of the Lorent z factor along 
a force-free field line ([Mest el fc Shibatalll994l : lOontopoulosI 
los et al.l 



1995; ContoDoulos et al. 1999). For an stationary, axisym- 
metric model the energy (Bernoulli) equation alone deter- 
mines the field- aligned velocity and the frozen-in conditions, 



B r 



B° 



FU FL 



B4> 



apply to the fluid velocity and by uf*' B a = 0, 



B r 



fUem 



B* 



(45) 



(46) 
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apply to the electromagnetic velocity. These equations can 
be solved for the electromagnetic 3- velocity to obtain 



V E M 



B l B, 



( 'f + 



B 2 



(47) 



where and </> M are the time and cf> Killing vectors and _B M = 
g^vB" . That is, for stationary, axisymmetric solutions the 
magnetic field and Qf alone determine the electromagnetic 
velocity, unlike in general as given by equation 118H . 



3 ALGORITHM TESTS 

A parameter space test of the analytic GRFFE inversion 
described above was performed to evaluate the range of al- 
lowed Lorentz factors the inversion is capable of handling. 
The procedure is to start with a known v % and B 1 for a 
specific point in space-time, determine T^, and then invert 
to get v l . Double precision is used for all quantities. The 
typical failure mode is that the velocity is determined to be 
space- like, and this is nearly coincident with a significant 
increase in the error in the inversion. Of interest is the gen- 
eral maximum value of the Lorentz factor below which no 
failures occur. 

In Minkowski space-time tests, the maximum Lorentz 
factor is u* ~ 10 10 before roundoff error causes the inversion 
to suddenly have significant error. Below this it', the rela- 
tive error is similar to machine error. As an extreme test, the 
space-time point is chosen in Kerr-Schild coordinates with 
a Kerr spin parameter of a/M — 0.9375 for a point on the 
horizon at 6 — 7r/4. For grid-aligned flows the maximum 
Lorentz factor is u l ~ 10 10 as before. However, for arbi- 
trary flow directions, the inversion can fail for u l > 2200. 
Most astrophysical flows of interest have u < 2000, so this 
is sufficient for our purposes. Otherwise a smaller machine 
precision should be used. 

The GRFFE formulation is coupled to HARM to test 
the formulation's ability to handle typical force-free prob- 
lems. The GRFFE formulation is used in HARM to per- 
form the Minkowski space-time test ca lculations that are de- 
scribed bv iKomissarovl J2002bl l2004a|) . These tests include: 
1) a fast wave ; 2) a degenerate Alfven wave ; 3) a three- 
wave problem ; 4) a problem that evolves to B 2 — E 2 < ; 
5) a standing Alfven wave ; and 6) two current sheet prob- 
lems. Since our GRFFE formulation assumes E-B = 0, their 
non-degenerate Alfven wave test is not considered. 

Their notation of the "wave frame" fields E' and B' 
are equivalent to the HARM notation for = and & M for 
E' — 0, and otherwise the wave frame, E' , and B' can be 
processed through the inversion routine to setup an initial 
lab-frame B 1 and v l from their E' , B' , and speed (i. The 
other problems they setup only require the lab-frame E and 
B, and again our GRFFE inversion can be used to obtain 
the initial lab-frame B % and v t as long as the problem is 
degenerate. Note that our sign conventions for E and B 
agree with theirs. We estimate that they used about 200 
grid zones for their tests and so use the same number for 
our tests. We use the same final time, box size, and plot 
labels for easy comparisons. W e shift the x = position t o 
be similar for each set of tests in lKomissarovl (l2002bl . r2004af) . 
A Courant factor of 0.9 and the monotonized central limiter 



is used for all tests except an additional current sheet test 
with parabolic spatial interpolation. 

Notice t hat apart fro m the initial condition given in 
IKomissarovl J2002bt [2004ah . for the fast wave one requires a 
relation determined earlier in their paper: E z = C — i-ifB v , 
where C is constrained such that B 2 — E 2 < (C = 1 was 
chosen) and (if — +1 was chosen by them. 

Fi gure ^ show s the suite of tests demonstrated in 
IKomissarovl (120 02b ) . For the top 3 panels, the solid line de- 
notes the analytic solution at the initial and final time. The 
diamonds denote the numerical solution at the final time. 
The top panel shows B2 for t he fast wave. The he ad of the 
wave is more resolved than in IKomissarovl i2002rj) . but the 
tail of the wave is slightly less resolved. The second and 
third panels from the top show B2 and B3 for the degen- 
erate Alfven wave problem. The code does well to capture 
both fast and Alfven waves. The next two panels show B2 
and Ei for th e three-wave problem . The waves are resolved 
similarly as in lKomissarov ( 2002b). The bottom panel shows 
b 2 /B 2 = (B 2 - E 2 )/B 2 for the smoothed B 2 - E 2 -> test 
problem. Overall, the GRF FE version o f HARM performs 
comparably to the code by IKomissarovl j2002bl) based on 
more complicated Riemann solvers. 

F igure [5] shows the suite of tests demonstrated in 
IKomissarovl (2004a). The top panel shows B z for a station- 
ary Alfven wave. T his wave is much more resolved than in 
IKomissarovl (|2004a) and indicates that the code's effective 
resistive diffusion coefficient is quite low. The next panel 
sh ows E z and B v for a current sheet problem as described 
in lKomissarovl j2004al) with Bo = 0.5. The features are well- 
resolved. The two bottom panels show the second current 
sheet problem with Bq — 2. The upper of the two is using the 
MC limiter, w hile the lower of the two is using a parabolic 
interpolation (IColelralll98l . which gives similar results to 
IKomissarovl i2004aF "ln the fast wave region the MC lim- 
iter gives a Lorentz factor of T = 1.90, while the parabolic 
method gives F = 2000, the largest allowed Lorentz factor. 
For either method, the region at x 1 — 0.5 within the current 
sheet reaches B 2 — E 2 ~ 0, but the Lorentz factor limiter 
(here limited to ymax = 2000) keeps the code stable despite 
the presence of the current sheet. No complicated dissipation 
model had to be included to achieve such a result. 



4 PHYSICAL MODELS 

This section considers astrophysical models for which the 
GRFFE approximation is a reasonable one. The GRFFE 
formulation is used in HARM. 



IKomissarovl <2004al) study the Wald ( Waldlll974h and 
split-monopole BZ77 jBlandford fc Znaiek ^)77T) solutions 



for slowly and moderately rapidly rotating black holes. Of 
particular interest is whether the Wald and split-monopole 
solutions can be bet ter represent e d comp ared to as shown 
in figures 3 and 4 of Komissarov ( 2004a|) . We consider the 
monopole and split-monopole for slowly rotating black holes 
in this paper. 

The other m odels are considered in a separate paper 
iMcKinnevl2005efl . which demonstrates our GRFFE formu- 
lation's ability to handle the current sheet in the a ctual 
spto-monopole solution of iBlandford fc Znaiekl il977f) and 
the Wald problem with a rapidly rotating black hole spin. In 
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Figure 1. Suite of tests described bv IKomissarovl |2002b). Solid curves represent the initial and final analytic solution for the top three 
panels. Diamonds represent the numerical solution at the final time, except for the plot of b 2 that shows both times. Tests are as follows 
from top to bottom: 1) B2 for fast wave ; 2) B2 and B3 for degenerate Alfven wave ; 3) B2 and E\ for three-wave problem ; and 4) 
initial and final smoothed b 2 = B 2 — E 2 — » problem. 



that paper fo r the Wald p r oblem^ we were able to reach sim- 
ilar solutions iKomissarovl i2005al) who used an MHD model 
to avoid significant reconnection in the current sheet t hat 
developed in their force- free models Komissarov ( 2004a) . 

Neutro n star magnetosp heres are studied in a sepa- 
rate paper iMcKinn ev 2005d), which demonstrates that the 
GRFFE code can be used to study pulsar magnetospheres 
even in the presence of a current sheet. That pa per shows we 
are abl e to avoid the problems encountered bv IKomissarovl 
(2005b) with force- free and the current sheet. We found sim- 
ilar results to, but more accurate than, their ideal MHD 
results. 



In this paper we focus on slowly rotating black hole 
magnetospheres that require general relativity and so 
full GRFFE. The interesting quantities that are plotted 
throughout the following sections are Qf = Fte/Fg,/,, which 
is also Qp — Ftr/F rt f, for a stationary, axisymmetric flow. 
This quantity often appears as a ratio to the black hole an- 
gular velocity of Qh = a/(2r+). The radial an d magnetic 
field strengths for these plots is defined as in IKomissarovl 
J2004al) . with B 1 = F . Also interesting is the conserved 
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Figure 2. Suite of tests described by Komissarov (20043). Solid curves represent the initial and final analytic solution except in the 
bottom panel. Pluses and diamonds represent the numerical solution at the final time. Tests are as follows from top to bottom: f) B z 
for stationary Alfven wave ; 2) E z and B y for current sheet with Bo = 0.5 ; 3,4) E z , b 2 , and B y for current sheet with Bo = 2. Third is 
with 2nd order method and fourth is with 4th order method. 



toroidal flux of sT^gBj, = = *F t<t> = y/=g~F re } This 
is because the electromagnetic energy flux is F% = — Tt = 
—B r Q.FB<f, and the electromagnetic angular momentum flux 
is Fl = Fe/Qf. Also of interest is the magnetic vector po- 
tential (A$), whose contours are plotted and represent flow 
surfaces. For axisymmetric, stationary flows these surfaces 
define surfaces of constant Qp, B^, F^/B 1 , and Fl/B 1 in 
dissipationless force-free electrodynamics. Finally, also of in- 
terest are the light surfaces, which are defined for the case 

1 Notice the slight change, for simplicity, in notation for B^ and 
B % from this point onward. 



of a purely rotational velocity fi = Qf leads to a null tra- 
jectory, i.e. from u a u a = —1, 

g tt + 2(lg^ t + Q?g H = 0. (48) 

4.1 Black Hole Magnetosphere: Blandford-Znajek 
Monopole and Split-Monopole 

iKomissarovl (|2002b) demonstr ated the stability of the pure 
monopole solution of iBlandford fc Znaiekl ll977l) by study- 
ing one hemisphere. Our GR FFE formulation ge n erates 
quantitatively similar results. iKrasnopolskv et alJ (120051) 
also use a method based on HARM and study the depen- 
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dence of the field (purely monopolar-type) geometry and 
field rotation frequency on black hole spin. They also stud- 
ied how the energy output varies with black hole spin up 
to a/M = 0.9999. They find that the energy output fol- 
lows E oc Qjj up to a/M = 0.98. This is in contrast to 
the energy output in the presence of a thick disk that fol- 
lows E oc Q% oc fl 3 p wh en accounting for the mass accretion 
rate iMcKinnevll2005aD . This suggests that the presence of 
matter can be important. However, there may be systems 
that are d ominated by a magnetosphere rather than a dis k 
(see, e.g.. Ilgumenshchev et alJl2003t IrTaravan et alJl20ojl . 
For such systems, the force-free limit may be sufficient to 
describe the energetics and geometry of the field lines. 

For the purposes of testing, the pure monopole 
Bl andford-Znaj e k solut ion is considered, as in section 5.2 
of iKomissarovl i2004al) . An identical numerical setup to 
lKorruss£a"ovl ~ il2004ar) is used. A grid is chosen with an 
outer radius of r = 260GA//c 2 . There are 150 equally- 
logarithmically-s paced radial z ones an d 100 uniform 8 zones. 
As mentioned in IKomissarovl i2004al) . the final result is in- 
sensitive to the details of the initial conditions, so the simple 
a/M = monopole solution is chosen as an initial condition. 
The evolution proceeds until t = 50GM/c 3 . 

Figure|21shows Qf/Qh (top panel), — B$, (middle panel; 
pluses), analytic Blandford-Znajek (BZ) solution for — B^ 
(middle panel; solid line), — _B r /10 (middle panel; short 
dashed line), B 9 (middle panel; long dashed line), and E ■ B 
for 8 — n — 0.1 a nd 8 = 0.1. The figu re is directly compara- 
ble to figure 2 in IKomissarovl j2004al) . The value of Qf/^h 
only varies monotonically fr om 0.5013 near the equator to 
0.5009 at the poles. Unlike in IKomissarovl <l2004al) . the value 
does not rise again near the poles, which suggests HARM is 
properly resolving the coordinate singularity. Their value is 
up to 0.503 near the equator, which suggest HARM is doing 
marginally better. The middle panels shows that the code is 
very accurately reproducing the BZ solution. The solid line 
indicating the analytic solution is not visibly deviating from 
the crosses that mark the numerical solution. The bottom 
panel just advertises the fact that HARM is uniformly pre- 
serving E ■ B = to machine precision, where figure 2 in 
IKomissarov l (2004a) shows nonuniform errors of order 10 

IKomissarovl (|2004a) study the Wald and split-monopole 
BZ77 solution for slowly and moderately rapidly rotating 
black holes. Of particular interest is whether the Wald and 
split-monopole solutions can be bett er represented com - 
pared to as shown in figures 3 and 4 of lKomissarovl (l2004al) . 
This paper demonstrates our GRFFE formulation's ability 
to handl e the current sheet in the ac tual spte-monopole so- 
lution of iBlandford fc Znaiekl ll 19771). 

We use an identical setup to lKomissarovl j2004al) for the 
split-monopole case with a black hole with spin a — 0.1 2 . 
There are 100 uniform 8 zones with 80 radial zones such 
that dr/r = Const.. The model is axisymmetric and the 
inner boundary is placed inside the horizon and the outer 
boundary is at r = 29GM/c 2 . T he model is run till t = 
5GM/c 3 as in IKomissarovl (l2004al) . 



2 IKomissarovl {2004a) pointed out in their section 5.3 that a = 
0.1, but then their figure 3 caption says "Schwarzschild" . The 
a = and a = 0.1 models show the same results, so this issue is 
not crucial. 



Figure 0] shows a two panel plot th a t can b e directly 
compared with figure 3 in IKomissarovl J2004a). In their 
model the current sheet rapidly reconnects even over this 
very short time period. Our results show essentially zero re- 
connection. This is primarily due to the modifications of the 
drift velocity described in section 12.51 and secondarily due 
to the low diffusivity of the algorithm as demonstrated in 
section [3] Without the drift velocity modificati on, however, 
our res ults are similar to shown in figure 3 of IKomissarovl 
J2004al) . So this demonstrates the usefulness of the relatively 
ad hoc approach of treating current sheets. 

For this model, we used the local Lax-Friedrich 
(LLAXF), rather than the Harten-Lax-van Leer (HLL), ap- 
proxim ate nonlinear Riemann solver (see, e.g. lGammie et all 
2003a). The HLL solver generated fluctuations at the inner- 
radial boundary that back propagate through the horizon 
due to numerical diffusion. HLL is explicitly causal. How- 
ever, in HARM and all numerical meth ods the auth or is 
aware of (for an exception see, e.g. lHawke e t al. 2005J), the 
stencil used to reconstruct the quantities to the Riemann in- 
terface is still acausal. For example, in HARM, the stencil is 
always centered and so connects information across the hori- 
zon depending upon the smoothness of the flow. The author 
knows of no numerical method that accounts for the char- 
acteristic inf ormation in constra ining the reconstruction to 
be causal. In iMcKinnevI <l2005el) . we consider modifications 
to the stencil that enforce strict causality. This eliminates 
all the problems described above. 



5 CONCLUSIONS 

We described a GRFFE formulation that allows GRMHD 
codes to directly evolve the GRFFE equations of motion. 
Rather than evolving the electric and magnetic fields, the 
velocity and magnetic field are directly evolved. This for- 
mulation strictly enforces E • B = 0, V • B = 0, and energy 
conservation, unlike prior codes. Established, accurate, and 
well-tested GRMHD codes can simply add a new inversion 
piece of code to their existing code, while continuing to use 
all the already-developed facilities present. 

We also introduced a simplified general model of the 
dissipation of the electric field to enforce the B 2 — E 2 > 
constraint. This limits the code to a regime of Lorentz fac- 
tors that the code can handle without significant numerical 
errors. 

A simplified general model was introduced to allow cur- 
rent sheets to be resolved, without reconnection, over many 
dynamical times. The other improvements to the code, such 
as strict enforcement ofE-B = 0, VB = 0, energy conserva- 
tion, and B 2 — E 2 > do not play a role in this improvement. 
Limiting the numerically induced drift velocity perpendicu- 
lar to the sheet was a crucial step to make such a force-free 
code use ful, and resolves the difficul ties encountered by prior 
authors llKomissarovll2004al . l2005bl) . For highly magnetized 
systems, our GRFFE results are as good as ideal MHD re- 
sults without the need to introduce an artificial evolution of 
the rest-mass and internal en ergy densities. Evi dence of this 
fact is shown in other papers (Mc Kinnevl2005c]ll3) . where we 
model particular astrophysical systems. 

Numerical tests showed that the GRFFE formulation 
as used in HARM is robust and accurate. HARM uses 
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Figure 3. Pure monopole Blandford-Znajek solution for black hole spin a/M = 0.1. Upper panel: Qp/flu- Middle panel: — (dotted: 
numerical, solid : analytic), B r (sho rt dashed), and B® (long dashed). Lower panel: E • B. Model at t = 50GM/c 3 . Directly comparable 
with figure 2 in IKomissarovl l2004aT) . 



a simplified nonlinear ap proximate Riemann solver as in 
iKrasnopolskv et alJ J2005I) . which, whi le being muc h simple r 
than the exact Riemann solution bv iKomissarovl |2002b) , 
produced as accurate results. 

For the pure monopole Bla ndford-Zna.je k mode l, 
we found similar resu lts as IKomissarovl ll2002bl) : 
IKrasnopolskv et alJ J2005I) . We also demonstrated our 



code's ability to handle cu r rent sh eets, which w as found to 
be diffi cult in IKomissarovl j2004al) . Here and in IMcKinnevI 
the split-monopole problem was able to be solved 
without significant reconnection. This allows one to study 
general magnetospheres with arbitrary cu r rents i n a sheet, 
such as the paraboloidal field. IMcKinnevI l)2005eh also dis- 
cusses the Wald solution, and it shows that our code is able 
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Figure 4. Contours of the magnetic vector potential component giving the field lines for the split-monopole Blandford-Znajck 
soluti on for black hole spin a/M = 0.1. Left panel: Initial Model. Right panel: Model at t = 5GM/c 3 . Directly comparable with figure 3 
in iKomissarovl l2004aT) . 



to obta in simi lar results as the M HD code of IKomissarovl 
J2005al) . while IKomissarovl <2004al) encountered difficulties 
with the current sheet reconnecting too fast. 
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